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Abstract 

A multigrid method for the Stokes system discretized with an ^"-conforming 
discontinuous Galerkin method is presented. It acts on the combined velocity and 
pressure spaces and thus does not need a Schur complement approximation. The 
smoothers used are of overlapping Schwarz type and employ a local Helmholtz 
decomposition. Additionally, we use the fact that the discretization provides nested 
divergence free subspaces. We present the convergence analysis and numerical 
evidence that convergence rates are not only independent of mesh size, but also 
reasonably small. 


1 Introduction 

The efficient solution of the Stokes equations is an important step in the develop¬ 
ment of fast flow solvers. In this paper we present analysis and numerical results 
for a multigrid method with subspace correction smoother, which performs very 
efficiently on divergence-conforming discretizations with interior penalty. We ob¬ 
tain convergence rates for the Stokes problem which are comparable to those for 
the Laplacian. 

Multigrid methods are known to be the most efficient preconditioners and 
solvers for diffusion problems. Nevertheless, for Stokes equations, the divergence 
constraint makes the solution process more complicated. A typical solution em¬ 
ploys the use of block preconditioners, e. g. 1 1 31122II27II281 . but their disadvantage 
is, that their performance is limited by the inf-sup constant of the problem. This 
could be avoided, if the multigrid method operated on the divergence free sub¬ 
space directly, and thus would not have to deal with the saddle point problem at 
all. Such methods have been developed in different context and have proven very 
successful as reported for instance by Hiptmair m for Maxwell equations and by 
Schoberl |33l for incompressible elasticity with reduced integration. 
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The main ingredients into such a method are a smoother which operates on 
the divergence free subspace and a grid transfer operator from coarse to fine mesh 
which maps the coarse divergence free subspace into the fine one. The second 
objective can be achieved by using a mixed finite element discretization for which 
the weakly divergence free functions are point-wise divergence free. For such a 
discretization, the natural finite element embedding operator from coarse to fine 
mesh does not increase the divergence of a function. Discretizations of this type 
are available, such as for instance in Scott and Vogelius |36ll39l . Neilan and coau¬ 
thors I I 141 161 and Zhang 1411421 . Here, we focus on the divergence conforming dis¬ 
continuous Galerkin (DivDG) method of Cockburn, Kanschat, and Schotzau tm 
due to its simplicity. 

Following the approach by Schoberl [33], in order to study smoothers for the 
Stokes equations, we first consider a problem on the velocity space only with 
penalty for the divergence. This leads to a singularly perturbed problem with an 
operator with a large kernel. When it comes to smoothers for such operators, there 
are two basic options. One approach is to smooth the kernel space explicitly, as 
proposed for instance by Hiptmair m and Xu in ED. The other option was pre¬ 
sented by Arnold, Falk, and Winther in 12151 and smoothens the kernel implicitly, 
while never employing an explicit description of it. 

We follow the implicit approach and use the same domain decomposition prin¬ 
ciple (i.e additive and multiplicative Schwarz methods and vertex patches), but 
instead of the Maxwell or divergence dominated mass matrix as in 12131 apply it to 
the DivDG Stokes discretization. Then, we prove the convergence of the multigrid 
method with variable V-cycle algorithm for the singularly perturbed problem. The 
second pillar we rest on is the equivalence between singularly perturbed, diver¬ 
gence dominated elliptic forms and mixed formulations established by Schoberl 
in Ii 32ll33| . This equivalence allows us to apply the smoother to a mixed formula¬ 
tion of nearly incompressible elasticity and then to proceed to the Stokes limit. 
As far as we know, the combination of these techniques has not been applied 
the DivDG method in m Since our analysis is based on domain decompo¬ 
sition, fundamental results are also drawn from the seminal paper by Feng and 
Karakashian m on domain decomposition for discontinuous Galerkin methods 
for elliptic problems. 

There is a close relation between our technique and the smoother suggested 
by Vanka in [38] for the MAC scheme: the MAC scheme can be considered the 
lowest order case of the DivDG methods (see |[23)). In this case, the subspace 
decomposition structure of Vanka smoother corresponds to Neumann problems 
on cells, while our smoother is based on Dirichlet problems for vertex patches. 
Generalizations of the Vanka smoother have been applied successfully to different 
other discretizations albeit their velocity-pressure spaces are not matched in the 
sense of □ (see for instance |37||40| and literature cited there). 

Recently, an alternative preconditioning method for Stokes discretizations of 
the same type as here has been introduced in {5] by Ayuso et al. Their method 
is based on auxiliary spaces introduced by Hiptmair and Xu in ED. The exact 
sequence property of the divergence-conforming velocity element plays a crucial 
role as in our scheme, but their preconditioner uses a multigrid method for the 
biharmonic problem to solve the Stokes problem. As a consequence, it is not 
possible to use the preconditioning method for no-slip boundary conditions. On 
the other hand, it has been demonstrated in 126] that the multigrid method here can 
be lifted to the biharmonic problem, providing an efficient method for clamped 
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boundary conditions. 

The paper is organized as follows. In Section [2] we present the model prob¬ 
lem and the DG discretization. The multigrid method and domain decomposition 
smoother are derived in Section[3] Section[4]is devoted to the convergence analysis 
of our preconditioning technique with the man result in Theorem[]]on page[9] The 
paper concludes with numerical experiments in Section [5] 


2 The Stokes problem and its discretization 

We consider discretizations of the Stokes equations 

—Am + Vp = / in 0, 

V-u =0 in fl, (1) 

m = u B on <9f2, 

with no-slip boundary conditions on a bounded and convex domain Q. C R d 
with dimension d = 2, 3. The natural solution spaces for this problem are V = 
Hq(Q; R d ) for the velocity u and the space of mean value free square integrable 
functions Q = Lq(Q) for the pressure p, although we point out that other well- 
posed boundary conditions do not pose a problem. 

In order to obtain a finite element discretization, we partition the domain ft 
into a hierarchy of meshes {Tt}i=o,...,L of parallelogram and parallelepiped cells 
in two and three dimensions, respectively. In view of multilevel methods, the index 
l refers to the mesh level defined as follows: let a coarse mesh To be given. The 
mesh hierarchy is defined recursively, such that the cells of T^+i are obtained by 
splitting each cell of Ti into 2 d congruent children (refinement). These meshes are 
nested in the sense that every cell of Ti is equal to the union of its four children. 
We define the mesh size hi as the maximum of the diameters of the cells of T i. 
Due to the refinement process, we have hi = 2~ e ho. 

By construction, these meshes are conforming in the sense that every face of 
a cell is either at the boundary or a whole face of another cell; nevertheless, local 
refinement and hanging nodes do not pose a particular problem, since they can be 
treated following | ; 20]l2Tl l. By we denote the set of all faces of the mesh T i, 
which is composed of the set of interior faces FJ and the set of all boundary faces 
Ff. 

We introduce a short hand notation for integral forms on Ti and on F« by 

(0.VOt £ = Y / (c/),ip) Fe = Y / fields, 

TgT £ "' T FG Ff F 

IMIt 4 = (Z) / w 2dx ) - Mf £ = (IZ / M 2ds ) ’ 

'T GT^ T ' 'F£W e jF ' 

The point-wise multiplication operator <j> © t[> refers to the product <jnp, the scalar 
product (p • %j) and the double contraction <j> : %[> for scalar, vector and tensor argu¬ 
ments, respectively. The modulus \<j>\ = \/(j> Q <f is defined accordingly. 

In order to discretize CD on the mesh T«, we choose discrete subspaces Xt = 
Vi x Qe, where Qi C Q. Following CD , we employ discrete subspaces Vi of the 
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space _ffo lv (f2), where 

H iw (n) = {»e L 2 (n-,R d )\v-v e l 2 (h)}, 

H^(Q) = {v £ H A "{Q)\vn = 0 on<9H}. 

Here, we choose the well-known Raviart-Thomas space |[30j, but we point out 
that any pair of velocity spaces Vt and pressure spaces Qt is admissible, if the key 
relation 


V-Ve = Qe ( 2 ) 

holds. The details of constructing the Raviart-Thomas space follow. 

Each cell T £ It can be obtained as the image of a linear mapping T't of 
the reference cell T = [0, l] d . On the reference cell, we define two polynomial 
spaces: first, Qk , the space of polynomials in d variables, such that the degree with 
respect to each variable does not exceed k. Second, we consider the vector valued 
space of Raviart-Thomas polynomials 14 = Qk + xQk■ Polynomial spaces Vt 
and Qt on the mesh cell T are obtained by the pull-back under the mapping 'Tt 
( see for instance El)- The polynomial degree k is arbitrary, but chosen uniformly 
on the whole mesh. Thus, we will omit the index k from now on. Concluding this 
construction, we obtain the finite element spaces 

V t = {v€ Jf 0 d,v (n)|VT 6 Te : v\ T £ V T }, 

Qt — {q £ Lq(Q)\\/T £ Tt : q\T £ Qt}- 


2.1 Discontinuous Galerkin discretization 


While the fact that Vt is a subspace of Hq"(Q) implies continuity of the normal 
component of its functions across interfaces between cells, this is not true for tan¬ 
gential components. Thus, Vt (jL H 1 (fi; R d ), and it cannot be used immediately to 
discretize 0. We follow the example in for instance mmm and apply a DG 
formulation to the discretization of the elliptic operator. Here, we focus on the in¬ 
terior penalty method 01291 . Let Ti and T2 be two mesh cells with a joint face F, 
and let it 1 and U 2 be the traces of a function u on F from Ti and T 2 , respectively. 
On this face F, we introduce the averaging operator 


„ T, tti + W2 
W 

In this notation, the interior penalty bilinear form reads 

a e {u, v) = (Vu, Vv) T( + 4(aigu 0 n}, fu ® nf } Fi 
- 2 (§Vi4, gn ® «}) F j - 2({Vu}, gn 0 

+ 2 (a L u,v) ¥ g - (d n u,v) ¥ g - (d n v,u) ¥ g. 


(3) 


(4) 


The operator “ 0 ” denotes the Kronecker product of two vectors. We note that the 
term 4{it 0 n J : -j[i) 0 n jj- actually denotes the product of the jumps of u and v. 

The discrete weak formulation of 0 reads now: find (ue,pt) £ Vt x Qt, such 
that for all test functions vt £ Vi and qt £ Qt there holds 


At 



at(ut,ve) + (pt,V-vt) - ( qt,V-ut) = F(vt,qt) = ( f,vt ). 

(5) 
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Discussion on the existence and uniqueness of such solutions can be found for 
instance in mmmm. Here, we summarize, that is symmetric. If <jl is suf¬ 
ficiently large, the form ae(., .) is positive definite independently of the multigrid 
level t £ [0, L], and that thus we can define a norm on Ve by 


v t \\ v = 


( 6 ) 


In order to obtain optimal convergence results and to satisfy Proposition l2.2l below. 
a l is chosen as a/hL, where h l is mesh size on the finest level L and <r is a 
positive constant depending on the polynomial degree. By this choice, the bilinear 
forms on lower levels are inherited from finer levels in the sense, that 

at{ue,ve) = a L (u e ,ve), Vui,vteVe. (7) 

A particular feature of this method is (see 11011111 1. that the solution ue is in 
the divergence free subspace 

Vt = {yi £ V t \V-vt = 0}, (8) 

where the divergence condition is to be understood in the strong sense. 

Proposition 0.1 (Inf-sup condition). For any pressure function q £ Qe, there 
exists a velocity function v £ Ve, satisfying 


inf sup 


(g, v-tt) 


q€Qe „£ y { k) 


> ye > 0 


(9) 


v e l 


' Qe 


where ye = c^J = cs/2 l ~ L and c is a constant independent of the multi grid 
level L 


Proof. The proof of this proposition can be found in [[35] Section 6.4], Indeed, a 
different result is proven there, with ye ~ 1 /k, where k is the polynomial de¬ 
gree in the /ip-method. Thorough study of the proof though reveals, that this 
fc-dependence is due to the penalty parameter of the form ae « k 2 /he. In our 
case, the penalty parameter depends on the fine mesh, not on he, such that ae ~ 
(he/h.L)/he, and that the role of the k 2 in the penalty is taken by the factor 
he/h L . □ 

For any u £ Ve, we consider the following discrete Helmholtz decomposition: 


U I _L 

u = u + u 


( 10 ) 


where u° £ V® is the divergence free part and u x belongs to its ae (., .)-orthogonal 
complement. For functions from this complement holds the estimate: 

Lemma 0.1. Let u x € Ve be ae(., .)-orthogonal to Ve, that is, 
ae{u^,v) = 0 V v £ Ve ■ 


Then, there is a constant a > 0 such that 
a 


d 2 


V-tr 


< ae{u 


m x ) < — 
h 


V-u 


( 11 ) 


ye is the inf-sup constant from inequality ©. 
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Proof. On the left side, we already argued above that ol is chosen large enough 
such that ae{.,.) is uniformly positive definite. Thus, we have with a positive 
constant a 

a l|Vw x || t 4 < a t ( u x ,u x ). 

But then, 

fv-u x , V-u xN ) < d 2 (Vm x , Vu 1 ) < —a^(ti x ,u x ), 

V J ct V Jt id 

On the right side, let q = V-ti x . Then q £ Qe due to 0. From ([9j, we conclude 
that there is u £ Ve such that V-w = q and ||n|| v < l/7£||g||. On the other 
hand, m x is the error of the orthogonal projection into V e °. Thus, m x must be the 
element with minimal norm, and in particular 11 u x 11 v < 11 m 11 v • □ 

2.2 The nearly incompressible problem 

We are going to prove convergence uniform with respect to the refinement level 
l of the proposed multigrid method for the Stokes problem by deviating twice. 
First, we provide estimates robust with respect to the parameter e of the nearly 
incompressible problem: find ( ue,pe ) £ Ve x Qe such that for all (vt,qe) £ 
Ve x Qe there holds 

At ((p*) > (9^)) + £ (P«’®) (12) 

This problem is connected with the simpler penalty bilinear form (see for in¬ 
stance also G3) 

Ae, s (u e ,ve) = a t {ue,ve) +e _1 (V-ue,S7-ve) (13) 

and the singularly perturbed, elliptic problem: find u,e £ Ve such that for all ve £ 
Ve there holds 

Ae,e{ue,ve) = (f,ve) ■ (14) 

Lemma 0.2. Let (M m ,p m ) be the solution to and u e be the solution to (H- 
Then, if 0 holds, the following equations hold true: 

Urn = Ue, and EPm = V -Urn =v -Ue- 

Proof. Testing CCD> with ve = 0 and qe £ Qe yields 

— (V-Um, qe) + e(pm,qe) = 0 Mqe£Qe- 

Due to 10- this translates to the point-wise equality ep m = V-u m . Substituting 
Pm in ( 112b and testing with the pair (ve, V -ve), which is possible again due to 0, 
we obtain that u m solves G3- 

If on the other hand u e solves 4H. we introduce p e = |V-u e , which translates 
to 

- (V-u e ,qe) + £(Pe,qe) = 0 V qe £ Qe, 

corresponding to d 1 2b tested with (0 ,qe). On the other hand, d 1 2b tested with 
(t7,0) is obtained directly from d 14b substituting p e . Thus, the equivalence is 
proven. □ 
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In order to help keeping the notation separate, we adopt the following conven¬ 
tion: the subscript e is dropped wherever possible. Furthermore, curly letters refer 
to the mixed form, while straight capitals refer to operators on the velocity space 
only. Thus: 

at (u, v) the vector valued interior penalty form 

Ae(u, v) the form of the singularly perturbed, elliptic problem <141 

^ the mixed bilinear form GU 

Similarly, capital letters like in Re for the smoother 1261 refer to the singularly 
perturbed, elliptic problem, while IZe is the corresponding symbol for the Stokes 
smoother 1241 . Additionally, we associate operators with bilinear forms using the 
same symbol: 

Ae, e '■ VeVe (Ae, s u, v) = At, e (u, v) = Ae(u, v) = Al{u, v) Vu,v€Ve 

At,e : Xt -» Xt (Ae, e x,y) = Ae,e {x,y) = Ae(x,y) = Ah(x,y) Vx,y G Xt 



3 Multigrid method 


In Section[2] we introduced hierarchies of meshes (Tf }. Due to the nestedness of 
mesh cells, the finite element spaces associated with these meshes are nested as 
well: 


Vo 

C 

Vi 

c . . 

. c 

Vl, 

Qo 

c 

Qi 

c . . 

. c 

Ql- 

Xo = Vo x Qo 

c 

X, 

c . . 

. c 

Vl x Ql = Xl 


This relation also extends to the divergence free subspaces, see for instance 1*261 : 

v 0 ° C V? C ... C V°. (15) 


The nestedness of the spaces implies that there is a sequence of natural injec¬ 
tions le : Xt —> Xe+i of the form It(ve, qe) = ( Ie, u ve, h, P qe), such that 

h,u ■ Ve —» Vt+i, It, v '■ Qe —> Qe+i, (16) 

h,u : Ve° -»■ Ve° +1 . (17) 

The L 2 -projection from A'^+i — > Xe is defined by qe) = {Ie,u v e> ^e, P 1e) 

with 

(yt+i - Ie, u ve+i,we) = 0 Mwe G Vt (qe+i - Ie, P qe+i,re) = 0 Mre G Qe- 

(18) 


The ,4-orthogonal projection Vt from (Vi x Ql) —> (Ve x Qe) is defined by 


Al [Ve 


= Al 


(19) 


for all ( u,p ) G (Vl x Ql), (ve,qe) G Ve x Qe- Similarly, The A-orthogonal 
projection Pe from Vl —¥ Ve is defined by 


A L (Peu,vt) = A L (u,ve) 


( 20 ) 


for all u £ Vl, ve £ Ve- 
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3.1 The V-cycle algorithm 

In this subsection we define V-cycle multigrid preconditioners Be, e and Be, s for 
the operators Ae,e and Ae lB , respectively. For simplicity of the presentation, we 
drop the index s. 

First, we define the action of the multigrid preconditioner Be : Xe — » Xe 
recursively as the multigrid V-cycle with m(£) > 1 pre- and post-smoothing steps. 
Let Tie be a suitable smoother. Let Bo = Aq 1 . For £ > 1, define the action of Be 
on a vector Ce = {ft, gt) by 

1. Pre-smoothing: begin with (uo,Po) = (0, 0) and let 


/ Ui 

VPi 


/ Mi-1 
\Pi-l 


+k ‘{ c ‘ ~ a ‘ fc:)) 


i = 1,..., m(£), (21a) 


2. Coarse grid correction: 


\Pm(£) + l J 


\Pm(£)) 


+ Be—lie—i 


^Ce — Ae 


\Pm(t)) 


(21b) 


3. Post-smoothing: 


/ Ui 


/ Mi-1 

\Pi-l 


+ Tie 


^Ce — At 



i — m(t) + 2,..., 2 m(£) + 1 
(21c) 


4. Assign: 


BeCt = 


( M2m(€) + l'\ 
\P2m(r) + l / 


(2 Id) 


We distinguish between the standard and variable V-cycle algorithms by the 
choice 

j m(L) standard V-cycle, 

ym(L)2 L ~ e variable V-cycle, 

where the number m(L) of smoothing steps on the finest level is a free parameter. 
We refer to Bl as the V-cycle preconditioner of Al- The iteration 

is the V-cycle iteration. 

The definition of the preconditioner Be : Ve —> Ve for the elliptic operator Ae 
follows the same concept, but dropping the pressure variables. 


3.2 Overlapping Schwarz smoothers 

In this subsection, we define a class of smoothing operators Tie based on a sub¬ 
space decomposition of the space Xe. Let Me be the set of vertices in the tri¬ 
angulation T^, and let Te lV be the set of cells in T t sharing the vertex v. They 
form a triangulation with N(N > 0) subdomains or patches which we denote by 
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The subspace Xt jtl = Vt, v x Qe tV consists of the functions in Xt with support 
in Qt tV . Note that this implies homogeneous slip boundary conditions on dQe lV for 
the velocity subspace Vt, v and zero mean value on Vlt >v for the pressure subspace 
Qe, v - The Ritz projection Ve,v '■ Xt —» Xt, v is defined by the equation 


A fe) 





£ Xl. v 


(23) 


Note that each cell belongs to not more than four (eight in 3D) patches Tt lV , one 
for each of its vertices. 

Then we define the additive Schwarz smoother 

7 Ze=r) ^2 'Pt, v AJ 1 (24) 


where r] £ (0,1] is a scaling factor, IZi is L 2 symmetric and positive definite. 

Similarly, we define smoothers of the singularly perturbed elliptic operator At, 
namely, Pi : Vi —> Vi rV is defined as 


At (Pt, v ut, vt, v ) = At ( ut , vt, v ) Vv t, v £ Vt, v , (25) 


and the additive Schwarz smoother is 


Rt — i 7 ^2 p t,vA t 1 ■ (26) 

4 Convergence analysis 

In this section, we provide a proof of the convergence for the variable V-cycle 
iteration with additive Schwarz preconditioning method. Our proof is based on 
the assumption that the domain O is bounded and convex, which will be omitted 
for simplicity in the statement of following theorems and propositions. Our main 
result is: 

Theorem 1. The multilevel iteration X — BlAl for the Stokes problem 0 with 
the variable V-cycle operator Bl defined in Section\Jd\employing the smoother 
IZi defined in equation ( 124b with suitably small scaling factor t] is a contraction 
with contraction number independent of the mesh level L. 

Proof. First, we consider the nearly incompressible problem G2}- For this weak 
formulation, we have by Theorem[3] that the multigrid method X — Bl,eAl,h is 
equivalent to the method I — Bl.eAl^ applied to the singularly perturbed prob¬ 
lem US in the velocity space. 

Convergence of the multilevel iteration / — Bl, s Al,e is shown in Theorem[2] 
for all e > 0 with a contraction number <5 < 1 independent of L and e. Thus, by 
Theorem [3] the same holds fori — Bl,eAl,e with positive e. 

Finally, in we can let e converge to zero. The limit yields the well-posed 
Stokes problem 0. and since the contraction number 5 is independent of e, we 
obtain uniform convergence with respect to the mesh level L in the limit e —» 
0 . □ 

The theorems and lemmas of the following subsections serve to establish the 
building blocks of the proof of Theorem[I| 
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4.1 The singularly perturbed problem 

Theorem 2. Let Re be the smoother defined in dj with suitably small scaling fac¬ 
tor r). Then, the multilevel iteration I — BlAl with the variable V-cycle operator 
Bl defined in Section \T7\ is a contraction with contraction number independent of 
the mesh level L and the parameter e. 

The proof of this theorem is postponed to page!12land relies on 
Proposition 2.1. If Re satisfies the conditions: 

Al ((/ — RtAe)w, w) > 0, \/w G Ve (27) 

and 

(Rf 1 [I - Pi- i]w, [/ - Pi- i]w) < PeA L {[I - Pi- i]u>, [/ - P t -i]w), Mw G V t 

(28) 

where fit = 0(3j) is defined in equation below. Then 

0 < Al((I — BeAe)w, w) < SAl(w, w), \/w£Ve (29) 

where S = and C are defined in Lemma \2B\ 

Proof. In the case of self-adjoint operators Ae which are inherited from a common 
bilinear form a(.,.), this proposition would be part of the standard multigrid theory 
if /3 1 were constant. Its proof can be adapted from similar theorems in 12H81191 We 
will prove the version needed here in the appendix. □ 


In the remainder of this section we use several propositions and lemmas to 
establish our smoother Ri satisfies the assumptions of Proposition 12. II For u G 
(I — Pi-i)w with arbitrary w G it follows from the discrete Helmholtz decom¬ 
position in Section[2]and the projection operator Pe, v in Section [L2l that u admits 
a local discrete Helmholtz decomposition 

Uv = uZ + Uv (30) 

Lemma 2.1. Given L 2 -symmetric positive definite Re defined in \3.2\ and sym¬ 
metric positive definite Al(-, ■) defined in , there exists a constant p G (0, 1] 

independent of I such that 


p(R e 1 u, u) 


inf 

u-u v 

Yl-uUv =u 


'y ^ Al(V'v, 

■v 


Uv') 


(31) 


Proof. The following proof can be found in 0 for the L 2 -inner product instead of 
ai(., .). We copy it here to ascertain that it does not depend on the actual structure 
of the operator Al since it is purely algebraic. Thus, it applies to the operator Al 
in this paper as it applies to the different operator applied there. Recall that 


Re = fl 5Z ]D LvAi 1 =ri ^2 Pi.vAl 1 - 


(32) 


From 


u = 


u v 

V (zNg 


(33) 
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we get 


r/(R e 1 w, u) = r\ Y^ (R e 1 u, u v ) 


= V Y, {AlPi^Aj^ 1 R t 1 u,u v ) 


(34) 

(35) 


vdzJfl 


<V 2 l E [AL^Pi^Ab 1 R e 1 u,A L 1 R l 1 u) > < Y, ( Alu v ,u v ) 


V(zzJ\fg 


I v€iJfg 


= tj 2 {(Aiti.A^ii/u )} 2 Y, (A l u v ,u v ) 


uGA/'^ 


(36) 

(37) 


= r / 2 {(«, R e 1 u)} 2 < Y, ( A L u v ,Uy ) > 

^ vGAfg J 

The above inequality works for arbitrary splitting, hence we have 

r){.RJ l u, u) < E Al(Uv,Uv) 

f ga/^ 

For the choice = Pt, v Pt Af 1 R~ 1 u we get 

p(Rf 1 u, u) = inf Y, A L (u v ,u v ) 
l V Ju. v tlv£*e 


(38) 


(39) 


(40) 


□ 


Lemma 2.2. Given the local Helmholtz decomposition in Go}. For any Uv € Ve, v 
, there exists constant C\ independent of multigrid level satisfying: 

E ll V ' w ^l| 2 < Ci E a e (Uv,Uv) (41) 

Proof. It follows from Lemma lOTl that the estimate || V■ u x || 2 < Cae(u ± ,u ± ) 
hold for all m x E Vi. It is easy to see that Ve, v is a subspace of Ve for any v, 
so the estimate are also valid on any patch. In 2-D case, one cell could at most 
be sharing by four patches(eight patches in 3D). Hence there exists a constant C\ 
independent of multigrid level such that the estimates holds for the summation of 
local estimates. □ 


Proposition 2.2. Given the overlapping subspace decomposition ofVi in \3.2\ and 
the interior penalty bilinear form ae(u, v) iit 0- Assume is chosen sufficiently 
large, the following estimate holds on each level i. Then, there is a constant C 2 
which is independent of multigrid level such that for any u £ Vi holds 

Yj cie(uv,u v ) < C 2 de(u, u) (42) 

■dGA/^ 
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Proof. For a fixed L, the penalty constant oi is a/h,L which is greater than 'a/he. 
For the latter, this is a standard result: the proof and details on the choice of a can 
be found in (Up. 1361], □ 

Proof of Theorem^ Recall the definition of Al -orthogonal projections Pi and 
Pi which restrict the projection on Qe lV ( zero elsewhere). Following [2j, we show 
that if 0 <? 7 < 1/4, the smoother Ri satisfies the first condition in Theorem[2 

For weVi 

Al{[I — ReAt]w,w) = Al{w,w) - r) Y. A l (Pi iV w,w) (43) 


but 

A L (Pe,vW,w) = Al(Pi, v w,Pi, v w) < Al(w, w) 5 Al(Pi, v w, Pe, v w)* (44) 
so 

Y. A l (Pi, v w,w) < Yj Al(w,w) < 4A l (w,w) (45) 

vga fg 

Flence the first hypothesis holds. 

Thus, it remains to check the second condition which could be reduced to the 
following problem: forw = (/ — Pt-\)w (where w £ Vf) with the decomposition 
u = ^ m„, there is a constant C such that 


Yj Al(u v ,u v ) < CA l (u, u) (46) 

v£ Mg 

Following Lemmas l0.1ll2.1ir2.2l and Proposition ^. 21 we get: 


^ ^ Al {uv . 

,u v ) = Y/ { a e( u v,Uv) + £ *(V • u v , V • «„)} 

(47) 





= Yj Ut>) +e _1 (V • u„, V • m^t)| 

(48) 


< C 2 at{u, u) + Ci—e^aiOuviUv) 

z ' a 

v£Afg 

(49) 


< C 2 at(u , u) + e~ 1 Ci—ae(u J ~ J u ± ) 
a 

(50) 


<C 2 aAu,u)+£~ 1 Ci-- (y-u ± ,S7-u ± ) 

«7 1 

(51) 


= C 2 aAu, u) + e^Ci — — (V-m, V-w) 

(52) 


< max <C 2 , Ci — — \ Al{u,u) 

l «7/J 

(53) 


= CeA L (u,u) 

(54) 

where Ci = max <j 


□ 
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Now set 


Pi = -Ci 

V 


(55) 


We have verified the two conditions in Proposition 12.11 

Lemma 2.3. Given Pi above and m(l) defined in 13.71 there is a constant C such 
that 


Pi 

2 m(£) 


< C 


(56) 


Proof. We will discuss this inequality in two cases: first, if Pi = then 


Pi ^ 

2 m(e) 2 m 0 2 L ~ e ~ 2m 0 


(57) 


On the other hand, if p = -C\ — ~ 

7 ' 7] cx 'y 


Pi 

2 m(t) 


Cl C _ 

2m 0 2 L ~ e 


Ci—~ Ciii 

_ q c ^ a c 

2m 0 i/2 1 - < ““ 2m 0 


=: C 


(58) 

□ 


4.2 The mixed problem 

Secondly, we will discuss the Stokes equation in mixed variables. Set Xi :S : = 
{( ui,pi ) £ Xi : V -ui = epi}. Now, it remains to show the equivalence between 
the multigrid algorithms. 

Proposition 2.3. The multigrid components fulfill the following properties: 

1. The smoother 7 Zi for the mixed problem defined in preserves Xi i6 . On 
the subspace it is equivalent to the smoother Re in primal variables. This 
means for ( ui,pi ) £ Xi iB and 

GO-*(2) 

there holds ( Ui,pe ) £ Xi lC and 

ue = Reue (60) 

2. The prolongation maps Xi- i iE into Xi tS . On the subspace it is equiva¬ 
lent to the prolongation Ie in primal variables. This means for (ui~ i ,pe~ i ) £ 
X(- i, e and 

(ui\ =x (ui. 1 \ 

\plj \pl-lj 

there holds ( ue,pi ) £ Xi iC and 

ue = Ie(ue-i) 


(61) 


(62) 
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3. The coarse grid solution operator maps Xp-i <E into Xe tS . On the subspace 
it is equivalent to the coarse grid solution operator in primal variables. This 
means for (ut,pt) £ Xe. e and 

, 63 ) 

there holds (ue~i, pe-i) £ Xe- i, e and 

ue-i — Aiut (64) 

Proof. The proof of this proposition can be found for the operators there in (33] p. 
93]. We do not provide it here since the arguments are purely linear algebra, and 
thus apply independent of the actual bilinear form. □ 

Theorem 3. The multigrid algorithm in mixed variables preserves the space X/ >E . 
On this subspace it is equivalent to the multigrid algorithm in primal variables. 
This means for ( ue,pe ) £ Xt, e and (ue,pe) = Be(ue,pe) there holds (ue,pe) £ 
Xe, E and 

ut = Beue (65) 

where Be and Be are the corresponding multigrid operators for each algorithm. 
Proof. The multigrid operator Be fulfills the recursion 

B 0 = Ao\ ( 66 ) 


Be = ( Re) me (I-Ie(I - (Be-^Af^Ie-if A e )(Re) me , (67) 

and the mixed operator Be fulfills a corresponding one. Then we apply the above 
proposition, and the theorem is proved by induction. □ 


5 Numerical results 

We test the additive Schwarz method which we have analyzed in the preceding 
section in order to ascertain that the contraction numbers are not only bounded 
away from one, but are actually small enough to make this method interesting. 
Furthermore, we conduct experiments, which go beyond our analysis, in particular 
regarding the choice of the penalty parameter and the number of smoothing steps 
on lower levels. 

The experimental setup for most of the tables is as follows: the domain is 
£7 = [— 1,1] 2 , the coarsest mesh To consists of a single cell T = Q. The mesh 
Te on level £ is obtained by dividing all cells in 'fe-i into four quadrilaterals by 
connecting the edge midpoints. Thus, a mesh on level l has 4 e cells, and the length 
of their edges is 2 1 ~ l . The right hand side is / = (1,1). For the relaxation 
parameter in the additive Schwarz method, we found that 0.5 is the value which 
provides the best results for all experimental setups, hence we keep it there in all 
the following experimental setups. 
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level 

RT\ 

rt 2 

3 

4 

4 

4 

4 

4 

5 

4 

4 

6 

4 

4 

7 

4 

4 

8 

4 

5 


Table 1: Number of iterations n§ to reduce the residual by 10 6 with the variable 
V-cycle algorithm with penalty parameter dependent of the finest level mesh size. 



m{£) 

= 1 

m(£) 

= 2 

level 

RTi 

rt 2 

RTi 

rt 2 

3 

7 

1 

4 

4 

4 

7 

7 

4 

4 

5 

7 

7 

4 

4 

6 

7 

7 

4 

4 

7 

8 

8 

4 

4 

8 

8 

8 

4 

4 


Table 2: Number of iterations n§ to reduce the residual by 10 6 with the standard 
V-cycle iteration with one and two pre- and post-smoothing steps. Penalty parameter 
dependent of the finest level mesh size. 


In Table Q] we first test the additive Schwarz smoother using variable V-cycle 
algorithm on a square domain with no-slip boundary condition. For the penalty 
constant in the DG form 0, we choose the penalty parameter as <j/Hl , where a = 
(fc + l)(fc + 2), on the finest level L and all lower levels l. Results for different pairs 
of RTk/Qk are reported in the table which show the fast and uniform convergence. 

In Table [2] we keep the same experimental setup and present iteration counts 
for the standard V-cycle algorithm with one and two pre- and post-smoothing steps, 
respectively. Although our analysis does not apply, we still observe uniform con¬ 
vergence results. We also see that the variable V-cycle with a single smoothing 
step on the finest level is as fast as the standard V-cycle with two smoothing steps, 
and thus the variable V-cycle is more efficient. 

In Table [3] we test the variable and standard V-cycles with penalty parameters 
depending on the mesh level £, namely a/he (where a is a positive constant de¬ 
pending on the polynomial degree) in the DG form 0. While our convergence 
analysis does not cover this case either, we observe convergence rates equal to the 
case with inherited forms. 

In Table 0 we provide results with GMRES solver and Bl as preconditioner 
for different experimental setups as in Tables Q] [2] and [^respectively. The second 
and third columns are results for variable V-cycle with penalty parameter depen¬ 
dent of the finest level mesh size. The fourth and fifth columns are the results for 
standard V-cycle with penalty parameter dependent of the finest level mesh size. 
The last two columns are the results for standard V-cycle with penalty parameter 
depend on the mesh size of each level. From this table, we see that the GMRES 
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variable 

standard 

level 


rt 2 


rt 2 

3 

4 

4 

7 

7 

4 

4 

4 

7 

7 

5 

4 

4 

7 

7 

6 

4 

4 

7 

7 

7 

4 

4 

7 

8 

8 

4 

5 

8 

8 


Table 3: Penalty parameter dependent on the mesh size of each level. Number of iter¬ 
ations n& to reduce the residual by 10 6 with variable and standard V-cycle iterations 
with m(L) = 1. 



variable 

standard 

noninherited 

level 

RT i 

rt 2 

RTi 

rt 2 

RT i 

rt 2 

3 

2 

2 

2 

2 

2 

2 

4 

3 

3 

3 

3 

3 

3 

5 

3 

3 

4 

3 

4 

4 

6 

3 

3 

5 

4 

5 

5 

7 

3 

3 

5 

5 

5 

5 

8 

5 

4 

6 

6 

8 

6 


Table 4: Number of iterations ne to reduce the residual by 1 0“ G with GMRES solver 
and preconditioner £>/,; variable and standard V-cycle with inherited forms, variable 
V-cycle with noninherited forms. One pre- and post-smoothing step on the finest level. 


method, as expected, is faster in every case. 


6 Conclusions 

In this paper, we have investigated smoothers based on the ones introduced by 
Arnold, Falk, and Winther for problems in H d "' in a variable V-cycle precondi¬ 
tioner for the Stokes system. We presented the convergence analysis and showed 
uniform contraction independent on the mesh level. In numerical experiments we 
showed that the contraction is not only uniform, but also very fast, thus making 
our method a feasible solver or preconditioner. 

In theory, the performance of the smoother relies on an exact sequence property 
of finite element spaces, in particular an i7 dlv -conforming discontinuous Galerkin 
discretization of the Stokes problem. Our experiments with the Taylor-Hood ele¬ 
ments, where the method fails, demonstrate that this is not an artifact of the analy¬ 
sis, but that the technique does not work due to the lack of an exact Hodge decom¬ 
position and nested divergence free subspaces. 
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A Proof for Proposition 12.11 

Following the proof in U, we want to show by induction on i that 

0 < A l ((I — BiAi)u,u) < SA l (u,u), \/u€Ve (68) 

For i = 1 is obvious since B\ = Af 1 . Now check if the above inequality hold 
for i = £ — 1. Recall the relaxation operator Ke = I — RiAe and the recurrence 
relation : 

I - BiAt = [(/ - P f _i) + (/ - B t -iAt-i)Pi-i] (69) 

The lower bound easily follows from the inductive hypothesis and the above iden¬ 
tity. For the upper bound, we use the induction hypothesis to obtain 


A l ((I - B e A e )u, u) < A l ([I - Pe-i\K™^u, K™ W u) + SA^Pj^K^u, K™ W u) 

(70) 

= (1 - S)A l ([I - Pe-^K^u, K™ w u) + 5A L (K™ W u, K™ w u). 

(71) 

Now by the orthogonality from OOt 


A l ({I - Pr-i]AT W M, [I - P«_i]A7 W u) (72) 

= A l ([I - P e -i]K™ W u,K™ W u) (73) 

= ([/ - P^i]K^ e) u, A t K™ W u) (74) 

= (R^ll - Pr_i]A7 W u, R e A e K™ m u) (75) 


< (R^ll - P t -i]K? (l) u, [I - Pe-i]K™ W u)? {R e A e K™ W u, A e K™ W u)i 

(76) 

< VFe{[I - Pe^K^u, [I - Pe-i]K™ W u)?(ReAeK™ W u, A e K™ W u)i 

(77) 

Hence, we get 

A l {[I - Pe-i]K™^u, I<™ W u) < /3t(RiAtK™ W u, A e K™ W u) (78) 

= 0iA l ([I - Ke-i\K 2 e mW u, u) (79) 

It follows from the positive semi-definiteness and d27l) that the spectrum of Ke 
is contained in the interval [0, 1]. Therefore, we have 

A l ([I- Ke-i]K 2rnW u,u) <A L {[I-K e -i]Kiu,u), for i<2m{£) 

(80) 

whence 

2m(£)-l 

A L ([7-^_i]A, 2ro(<) U ,w)< J2 Al{[I - Ke\Kiu,u) (81) 

= ^7Tv A ^ J - Ke]K 2 r W n, u) (82) 
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Combining CZQ) and m and following Lemma 12.31 we get 


A l ((I - B e At)u, u) < (1 - 5)^^A l ([I - K 2mW ]u , u) + 5A L (K™ W u, K™ W u) 

(83) 

< (l-5)CA L {[I - K 2rnW ]u,u) + 5A L {K™ w u,K™ w u) (84) 

= (1 - 5)CA l ( u,u) + [5 - (1 - 5)C}A l (K™ W u,K? W u) (85) 


The results now follows by choosing : 


S = (1-8)C, i. e„ 5 — ———— 
1 + C 


( 86 ) 
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